/* Copyright (c) Colorado School of Mines, 2011.*/
/* All rights reserved.                       */

/* SUKDSYN2D: $Revision: 1.11 $ ; $Date: 2011/11/12 00:40:42 $	*/

#include "su.h"
#include "segy.h"

/*********************** self documentation **********************/
char *sdoc[] = {
" SUKDSYN2D - Kirchhoff Depth SYNthesis of 2D seismic data from a	",
"             migrated seismic section					",
" 									",
"   sukdsyn2d  infile  outfile [parameters] 				",
"									",
" Required parameters: 							",
" infile=stdin		input migrated section				",
" outfile=stdout	file for output seismic traces  		",
" ttfile		file for input traveltime tables		",
"									",
" The following 9 parameters describe traveltime tables:		",
" fzt= 			first depth sample 				",
" nzt= 			number of depth samples 			",
" dzt=			depth interval 					",
" fxt=			first lateral sample 				",
" nxt=			number of lateral samples 			",
" dxt=			lateral interval 				",
" fs=			x-coordinate of first source			",
" ns=			number of sources				",
" ds=			x-coordinate increment of sources		",
"									",
" The following 6 parameters describe the migration section:		",
" fz=                   first z-coordinate in migrated section 		",
" dz=     		vertical spacing of migrated section 		",
" nz=           	number of depth points in migrated section	",
" fx=                   first x-coordinate of migrated section 		",
" dx=     		horizontal spacing of migrated section 		",
" nx=           	number of lateral points in migrated section  	",
"									",
" Optional Parameters:							",
" nt=501        	number of time samples				",
" dt=0.004      	time sampling interval (sec)			",
" ft=0.0        	first time (sec)				",
" nxo=1                 number of source-receiver offsets		",
" dxo=25                offset sampling interval  			",
" fxo=0.0               first offset  					",
" nxs=101               number of shotpoints  				",
" dxs=25                shotpoint sampling interval  			",
" fxs=0.0               first shotpoint 				",
" fmax=1/(4*dt)         maximum frequency in migration section (Hz)	",
" aperx=nxt*dxt/2  	modeling lateral aperature 			",
" angmax=60		modeling angle aperature from vertical		",
" v0=1500(m/s)		reference velocity value at surface		",
" dvz=0.0  		reference velocity vertical gradient		",
" ls=1	                flag for line source				",
" jpfile=stderr		job print file name 				",
" mtr=100  		print verbal information at every mtr traces	",
"									",
" Notes:								",
" This program takes a migrated seismic section and a set of travel time",
" tables generated using rayt2d for a specific background velocity model",
" and generates synthetic seismic data in the form of common shot gathers.",
" (Common offset gathers may be generated by using nxo=1.) (Demigration.)",
"									",
" This program is a tool which may be used for the migration residual	",
" statics estimation technique of Tjan, Audebert, and Larner 1994.	",
"									",
"1. The traveltime tables are generated by program rayt2d (or other ones)",
"   on relatively coarse grids, with dimension ns*nxt*nzt. In the	",
"   modeling process, traveltimes are interpolated into shot/geophone 	",
"   positions and migration section grids. 				",
"2. The input migration section must be an array of binary floats (no SU",
"   headers).								", 
"3. The synthesized traces are output in common-shot gathers in SU format.",
"4. The memory requirement for this program is about			",
"    	(ns*nxt*nzt+nx*nz+4*nr*nzt+3*nxt*nzt)*4 bytes 			",
"    where nr = 1+min(nxt*dxt,0.5*offmax+aperx)/dxo. 			",
"									",
NULL};
/*
 * Author:  CWP: Zhenyue Liu, 07/24/95,  Colorado School of Mines 
 *
 * References: 
 *
 * Tjan, T., F. Audebert, and K. Larner, 1994,
 *    Prestack migration for residual statics estimation in complex media
 *    (Appeared in 1994 Project Review, CWP-153.)
 *
 * Tjan, T., 1995, Residual statics estimation for data from structurally
 *    complex areas using prestack depth migration: M.Sc. thesis, Colorado
 *    School of Mines. (In progress.)
 *
 * Larner, K., and Tjan, T., 1995, Simultaneous statics and velocity
 *    estimation for data from structurally complex areas.
 *    (Appeared in 1995 Project Review, CWP-185.)
 *
 *
 * Trace header fields set: ns, dt, delrt, tracl, tracr, fldr, tracf
 *                          offset, sx, gx, trid, counit
 */

/**************** end self doc ***********************************/
 
/* prototypes of functions used internally */
void integ(float **mig,int nz,float dz,int nx,int m,float **migi);
void resit(int nx,float fx,float dx,int nz,int nr,float dr,
	float **tb,float **t,float x0);
void sum2(int nx,int nz,float a1,float a2,float **t1,float **t2,float **t);
void timeb(int nr,int nz,float dr,float dz,float fz,float a,
	float v0,float **t,float **p,float **sigb,float **cosb);
void maketrace(float *trace,int nt,float ft,float dt,
	float xs,float xg,float **mig,float **migi,float aperx,
  	int nx,float fx,float dx,int nz,float fz,float dz,
	int mzmax,int ls,float angmax,float v0,float fpeak,Wavelet *w,
	float **tb,float **pb,float **sigb,float **cosb,int nr,float **tsum,
	int nxt,float fxt,float dxt,int nzt,float fzt,float dzt);

/* parameters for half-derivative filter */
#define LHD 20
#define NHD 1+2*LHD

/* segy trace */
segy tr;

int
main (int argc, char **argv)
{
	int 	nt,nzt,nxt,nz,nxo,nxs,ns,nx,nr,is,ixo,ixs;
	int 	ls,jtr,mtr,tracl,mzmax;
	off_t nseek;
	float   ft,fzt,fxt,fz,fx,fxo,fs,fxs,dt,dzt,dxt,dz,dx,dxo,ds,dxs,
		ext,ezt,es,ex,ez,xo,s,xs,xg,fpeak;	
	float 	v0,dvz;
	float 	fmax,angmax,offmax,rmax,aperx;
	float **mig,**migi,***ttab,**tb,**pb,**cosb,**sigb,**tsum,**tt;
	
	char *infile="stdin",*outfile="stdout",*ttfile,*jpfile;
	FILE *infp,*outfp,*ttfp,*jpfp;
	Wavelet *w;


	/* hook up getpar to handle the parameters */
	initargs(argc, argv);
	requestdoc(1);

	/* open input and output files	*/
	if( !getparstring("infile",&infile)) {
		infp = stdin;
	} else  
		if ((infp=fopen(infile,"r"))==NULL)
			err("cannot open infile=%s\n",infile);
	if( !getparstring("outfile",&outfile)) {
		outfp = stdout;
	} else { 
		outfp = fopen(outfile,"w");
	}
	efseeko(infp,(off_t) 0,SEEK_CUR);
	warn("Got A");
	efseeko(outfp,(off_t) 0,SEEK_CUR);
	if( !getparstring("ttfile",&ttfile))
		err("must specify ttfile!\n");
	if ((ttfp=fopen(ttfile,"r"))==NULL)
		err("cannot open ttfile=%s\n",ttfile);
	if( !getparstring("jpfile",&jpfile)) {
		jpfp = stderr;
	} else  
		jpfp = fopen(jpfile,"w");

	/* get information for seismogram traces */
	if (!getparint("nt",&nt)) nt = 501;
	if (!getparfloat("dt",&dt)) dt = 0.004;
	if (!getparfloat("ft",&ft)) ft = 0.0;
	if (!getparfloat("fmax",&fmax)) fmax = 1.0/(4*dt);
	fpeak = 0.2/dt;
	if (!getparint("nxs",&nxs)) nxs = 101;
	if (!getparfloat("dxs",&dxs)) dxs = 15;
	if (!getparfloat("fxs",&fxs)) fxs = 0.0;
	if (!getparint("nxo",&nxo)) nxo = 1;
	if (!getparfloat("dxo",&dxo)) dxo = 50;
	if (!getparfloat("fxo",&fxo)) fxo = 0.0;
	
	/* get traveltime table parameters	*/
	if (!getparint("nxt",&nxt)) err("must specify nxt!\n");
	if (!getparfloat("fxt",&fxt)) err("must specify fxt!\n");
	if (!getparfloat("dxt",&dxt)) err("must specify dxt!\n");
	if (!getparint("nzt",&nzt)) err("must specify nzt!\n");
	if (!getparfloat("fzt",&fzt)) err("must specify fzt!\n");
	if (!getparfloat("dzt",&dzt)) err("must specify dzt!\n");
	if (!getparint("ns",&ns)) err("must specify ns!\n");
	if (!getparfloat("fs",&fs)) err("must specify fs!\n");
	if (!getparfloat("ds",&ds)) err("must specify ds!\n");
	ext = fxt+(nxt-1)*dxt;
	ezt = fzt+(nzt-1)*dzt;
	es = fs+(ns-1)*ds;

	/* check source and receiver coordinates */
 	for (ixs=0; ixs<nxs; ++ixs) {
		xs = fxs+ixs*dxs;
		for (ixo=0; ixo<nxo; ++ixo) {
			xg = xs+fxo+ixo*dxo;
			if (fs>xs || es<xs || fs>xg || es<xg) 
		err("shot or receiver lie outside of specified (x,z) grid\n");
		}
	} 

	/* get migration section parameters	*/
	if (!getparint("nx",&nx)) err("must specify nx!\n");
	if (!getparfloat("fx",&fx)) err("must specify fx!\n");
	if (!getparfloat("dx",&dx)) err("must specify dx!\n");
	if (!getparint("nz",&nz)) err("must specify nz!\n");
	if (!getparfloat("fz",&fz)) err("must specify fz!\n");
	if (!getparfloat("dz",&dz)) err("must specify dz!\n");
	ex = fx+(nx-1)*dx;
	ez = fz+(nz-1)*dz;
	if(fxt>fx || ext<ex || fzt>fz || ezt<ez) 
		err(" migration section is out of traveltime table!\n");

	if (!getparfloat("v0",&v0)) v0 = 1500;
	if (!getparfloat("dvz",&dvz)) dvz = 0;
	if (!getparfloat("angmax",&angmax)) angmax = 60.;
	if  (angmax<0.00001) err("angmax must be positive!\n");
	mzmax = dx*sin(angmax*PI/180.)/dz;
	if(mzmax<1) mzmax = 1;
	if (!getparfloat("aperx",&aperx)) aperx = 0.5*nxt*dxt;

	if (!getparint("ls",&ls))	ls = 1;
	if (!getparint("mtr",&mtr))	mtr = 100;

	fprintf(jpfp,"\n");
	fprintf(jpfp," Modeling parameters\n");
	fprintf(jpfp," ================\n");
	fprintf(jpfp," infile=%s \n",infile);
	fprintf(jpfp," outfile=%s \n",outfile);
	fprintf(jpfp," ttfile=%s \n",ttfile);
	fprintf(jpfp," \n");
	fprintf(jpfp," nzt=%d fzt=%g dzt=%g\n",nzt,fzt,dzt);
	fprintf(jpfp," nxt=%d fxt=%g dxt=%g\n",nxt,fxt,dxt);
 	fprintf(jpfp," ns=%d fs=%g ds=%g\n",ns,fs,ds);
	fprintf(jpfp," \n");
	fprintf(jpfp," nz=%d fz=%g dz=%g\n",nz,fz,dz);
	fprintf(jpfp," nx=%d fx=%g dx=%g\n",nx,fx,dx);
	fprintf(jpfp," \n");
	fprintf(jpfp," nxs=%d fxs=%g dxs=%g\n",nxs,fxs,dxs);
	fprintf(jpfp," nxo=%d fxo=%g dxo=%g\n",nxo,fxo,dxo);
	fprintf(jpfp," \n");
	
	/* compute reference traveltime and slowness  */
	offmax = MAX(ABS(fxo),ABS(fxo+(nxo-1)*dxo));
	rmax = MAX(es-fxt,ext-fs);
	rmax = MIN(rmax,0.5*offmax+aperx);
	nr = 2+(int)(rmax/dx);
	tb = ealloc2float(nzt,nr);
	pb = ealloc2float(nzt,nr);
	sigb = ealloc2float(nzt,nr);
	cosb = ealloc2float(nzt,nr);
	timeb(nr,nzt,dx,dzt,fzt,dvz,v0,tb,pb,sigb,cosb);

	fprintf(jpfp," nt=%d ft=%g dt=%g fpeak=%g \n",nt,ft,dt,fpeak);
	fprintf(jpfp," v0=%g dvz=%g \n",v0,dvz);
 	fprintf(jpfp," aperx=%g angmax=%g offmax=%g\n",aperx,angmax,offmax);
 	fprintf(jpfp," mtr=%d ls=%d\n",mtr,ls);
	fprintf(jpfp," ================\n");
	fflush(jpfp);

	/* allocate space */
	mig = ealloc2float(nz,nx);
	migi = ealloc2float(nz+2*mzmax,nx);
	ttab = ealloc3float(nzt,nxt,ns);
	tt = ealloc2float(nzt,nxt);
	tsum = ealloc2float(nzt,nxt);


	/* make wavelet */
	makericker(fpeak,dt,&w);

	/* set constant segy trace header parameters */
	memset((void *) &tr, 0, sizeof(segy));
	tr.trid = 1;
	tr.counit = 1;
	tr.ns = nt;
	tr.dt = 1.0e6*dt;
	tr.delrt = 1.0e3*ft;

	fprintf(jpfp," read traveltime tables \n");
	fflush(jpfp);

	/* compute traveltime residual	*/
	for(is=0; is<ns; ++is){
		nseek = (off_t) nxt*nzt*is;
		efseeko(ttfp,nseek*((off_t) sizeof(float)),SEEK_SET);
		fread(ttab[is][0],sizeof(float),nxt*nzt,ttfp);
		s = fs+is*ds;
		resit(nxt,fxt,dxt,nzt,nr,dx,tb,ttab[is],s);		
	}

	fprintf(jpfp," read migration section \n");
	fflush(jpfp);

	/* read migration section	*/
	fread(mig[0],sizeof(float),nx*nz,infp);

	/* integrate migration section for constructing anti-aliasing 
		filter	*/
	integ(mig,nz,dz,nx,mzmax,migi);
                       
	fprintf(jpfp," start synthesis ... \n");
	fprintf(jpfp," \n");
	fflush(jpfp);
	
	jtr = 0;
	/* loop over shots  */
	for (ixs=0,xs=fxs,tracl=0; ixs<nxs; ++ixs,xs+=dxs) {
		/* loop over offsets */
		for (ixo=0,xo=fxo; ixo<nxo; ++ixo,xo+=dxo) {
	    		float as,res;
	    		int is;
			xg = xs+xo; 
			/* set segy trace header parameters */
			tr.tracl = tr.tracr = ++tracl;
			tr.fldr = 1+ixs;
			tr.tracf = 1+ixo;
			tr.offset = NINT(xo);
			tr.sx = NINT(xs);
			tr.gx = NINT(xg);

	    		as = (xs-fs)/ds;
	    		is = (int)as;
			if(is==ns-1) is=ns-2;
			res = as-is;
			if(res<=0.01) res = 0.0;
			if(res>=0.99) res = 1.0;
			sum2(nxt,nzt,1-res,res,ttab[is],ttab[is+1],tsum);

	    		as = (xg-fs)/ds;
	    		is = (int)as;
			if(is==ns-1) is=ns-2;
			res = as-is;
			if(res<=0.01) res = 0.0;
			if(res>=0.99) res = 1.0;
			sum2(nxt,nzt,1-res,res,ttab[is],ttab[is+1],tt);
			sum2(nxt,nzt,1,1,tt,tsum,tsum);
				fflush(jpfp);

			/* make one trace */
			maketrace(tr.data,nt,ft,dt,xs,xg,mig,migi,aperx,
			  nx,fx,dx,nz,fz,dz,mzmax,ls,angmax,v0,fmax,w,
			  tb,pb,sigb,cosb,nr,tsum,nxt,fxt,dxt,nzt,fzt,dzt);
			
			/* write trace */
			fputtr(outfp,&tr);

	        	jtr++;
	        	if((jtr-1)%mtr ==0 ){
				fprintf(jpfp," generated trace %d\n",jtr);
				fflush(jpfp);
	    		}
		}
	}


	fprintf(jpfp," generated %d traces in total\n",jtr);


	fprintf(jpfp," \n");
	fprintf(jpfp," output done\n");
	fflush(jpfp);

	efclose(jpfp);
	efclose(outfp);

	    
	free2float(tsum);
	free2float(tt);
	free2float(pb);
	free2float(tb);
	free2float(cosb);
	free2float(sigb);
	free3float(ttab);
	free2float(mig);
	free2float(migi);

	return(CWP_Exit());
}


void integ(float **mig,int nz,float dz,int nx,int m,float **migi) 
/* integration of a two-dimensional array	 
  input: 
    mig[nx][nz]		two-dimensional array
  output:
    migi[nx][nz+2*m] 	integrated array 
*/
{
	int nfft, nw, ix, iz, iw;
	float  *amp, dw, *rt;
	complex *ct;


        /* Set up FFT parameters */
        nfft = npfaro(nz+m, 2 * (nz+m));
        if (nfft >= SU_NFLTS || nfft >= 720720)
                err("Padded nt=%d -- too big", nfft);

        nw = nfft/2 + 1;
	dw = 2.0*PI/(nfft*dz);

	amp = ealloc1float(nw);
	for(iw=1; iw<nw; ++iw) 
		amp[iw] = 0.5/(nfft*(1-cos(iw*dw*dz)));
	amp[0] = amp[1];

        /* Allocate fft arrays */
        rt   = ealloc1float(nfft);
        ct   = ealloc1complex(nw);

	for(ix=0; ix<nx; ++ix) {
        	memcpy(rt, mig[ix], nz*FSIZE);
       	 	memset((void *) (rt + nz), 0, (nfft-nz)*FSIZE); 
        	pfarc(1, nfft, rt, ct);

        	/* Integrate traces   */
		for(iw=0; iw<nw; ++iw){
			ct[iw].i = ct[iw].i*amp[iw];
			ct[iw].r = ct[iw].r*amp[iw];
		}

        	pfacr(-1, nfft, ct, rt);

        	for (iz=0; iz<m; ++iz)  migi[ix][iz] = rt[nfft-m+iz];
        	for (iz=0; iz<nz+m; ++iz)  migi[ix][iz+m] = rt[iz];
	}

	free1float(amp);
	free1float(rt);
	free1complex(ct);
}


/* residual traveltime calculation based  on reference   time	*/
  void resit(int nx,float fx,float dx,int nz,int nr,float dr,
		float **tb,float **t,float x0)
{
	int ix,iz,jr;
	float xi,ar,sr,sr0;

	for(ix=0; ix<nx; ++ix){
		xi = fx+ix*dx-x0;
		ar = abs(xi)/dr;
		jr = (int)ar;
		sr = ar-jr;
		sr0 = 1.0-sr;
		if(jr>nr-2) jr = nr-2;
		for(iz=0; iz<nz; ++iz)
			t[ix][iz] -= sr0*tb[jr][iz]+sr*tb[jr+1][iz];
	}
} 

/* sum of two tables	*/
  void sum2(int nx,int nz,float a1,float a2,float **t1,float **t2,float **t)
{
	int ix,iz;

	for(ix=0; ix<nx; ++ix) 
		for(iz=0; iz<nz; ++iz)
			t[ix][iz] = a1*t1[ix][iz]+a2*t2[ix][iz];
}
 
/* compute  reference traveltime and slowness	*/
      void timeb(int nr,int nz,float dr,float dz,float fz,float a,
	float v0,float **t,float **p,float **sig,float **cs)
{
	int  ir,iz;
	float r,z,v,rc,oa,rou,zc;


	if( a==0.0) {
		for(ir=0,r=0;ir<nr;++ir,r+=dr)
			for(iz=0,z=fz;iz<nz;++iz,z+=dz){
				rou = sqrt(r*r+z*z);
				t[ir][iz] = rou/v0;
				if(rou<dz) rou = dz;
				t[ir][iz] = rou/v0;
				p[ir][iz] = r/(rou*v0);
				sig[ir][iz] = v0*rou;
				cs[ir][iz] = z/rou;
			}
	} else {
		oa = 1.0/a; 	zc = v0*oa;
		for(ir=0,r=0;ir<nr;++ir,r+=dr)
			for(iz=0,z=fz+zc;iz<nz;++iz,z+=dz){
				rou = sqrt(r*r+z*z);
				v = v0+a*(z-zc);
				if(ir==0){ 
					t[ir][iz] = log(v/v0)*oa;
					p[ir][iz] = 0.0;
					sig[ir][iz] = 0.5*(z-zc)*(v0+v);
					cs[ir][iz] = 1.0;
				} else {
					rc = (r*r+z*z-zc*zc)/(2.0*r);
					rou = sqrt(zc*zc+rc*rc);
					t[ir][iz] = log((v*(rou+rc))
						/(v0*(rou+rc-r)))*oa;
					p[ir][iz] = sqrt(rou*rou-rc*rc)
						/(rou*v0);
				   	cs[ir][iz] = (rc-r)/rou; 
					sig[ir][iz] = a*rou*r;
				}
			}
	}
}


void maketrace(float *trace,int nt,float ft,float dt,
	float xs,float xg,float **mig,float **migi,float aperx,
  	int nx,float fx,float dx,int nz,float fz,float dz,
	int mzmax,int ls,float angmax,float v0,float fmax,Wavelet *w,
	float **tb,float **pb,float **sigb,float **cosb,int nr,float **tsum,
	int nxt,float fxt,float dxt,int nzt,float fzt,float dzt)
/*****************************************************************************
Make one synthetic seismogram 
******************************************************************************
Input:
**mig		migration section 
**migi		integrated migration section 
nt		number of time samples in seismic trace
ft		first time sample of seismic trace
dt		time sampleing interval in seismic trace
xs,xg		lateral coordinates of source and geophone 
aperx		lateral aperature in migration
nx,fx,dx,nz,fz,dz	dimension parameters of migration region
mzmax		number of depth samples in triangle filter
ls		=1 for line source; =0 for point source
w		wavelet to convolve with trace
angmax		migration angle aperature from vertical 	 
tb,pb,sigb,cosb		reference traveltime, lateral slowness, sigma 
		and cosine of emergent angle
nr		number of lateral samples in reference quantities
tsum		sum of residual traveltimes from shot and receiver
nxt,fxt,dxt,nzt,fzt,dzt		dimension parameters of traveltime table

Output:
trace		array[nt] containing synthetic seismogram
*****************************************************************************/
{
	int nxf,nxe,nxtf,nxte,ix,iz,iz0,izt0,nzp,jrs,jrg,jz,jt,jx,mz,iz1;
	float xm,x,dis,rxz,ar,srs,srg,srs0,srg0,sigp,z0,rdz,ampd,
	      sigs,sigg,coss,cosg,ax,ax0,pmin,
	      odt=1.0/dt,pd,az,sz,sz0,at,td,res,temp;
	float *zpt,**ampt,**ampti,**zmt,*amp,*ampi,*zm,*tzt,*work1;
	int lhd=LHD,nhd=NHD;
	static float hd[NHD];
	static int madehd=0;

	/* if half-derivative filter not yet made, make it */
	if (!madehd) {
		mkhdiff(dt,lhd,hd);
		madehd = 1;
	}
 
	/* zero trace */
	for (jt=0; jt<nt; ++jt)
		trace[jt] = 0.0;

	zmt = ealloc2float(nzt,nxt);
	ampt = ealloc2float(nzt,nxt);
	ampti = ealloc2float(nzt,nxt);
	amp = ealloc1float(nzt);
	ampi = ealloc1float(nzt);
	zm = ealloc1float(nzt);
	tzt = ealloc1float(nzt);
	zpt = ealloc1float(nxt);
	work1 = ealloc1float(nt);

	z0 = (fz-fzt)/dzt;
	pmin = 1.0/(2.0*dx*fmax);
	rdz = dz/dzt;

	xm = 0.5*(xs+xg);
	rxz = (angmax==90)?0.0:1.0/tan(angmax*PI/180.);
	nxtf = (xm-aperx-fxt)/dxt;
	if(nxtf<0) nxtf = 0;
	nxte = (xm+aperx-fxt)/dxt+1.0;
	if(nxte>=nxt) nxte = nxt-1;

	/* compute amplitudes 	*/
	for(ix=nxtf; ix<=nxte; ++ix){
		x = fxt+ix*dxt;
		dis = (xm>=x)?xm-x:x-xm;
		izt0 = ((dis-dxt)*rxz-fzt)/dzt-1;
		if(izt0<0) izt0 = 0;
		if(izt0>=nzt) izt0 = nzt-1;

		ar = (xs>=x)?(xs-x)/dx:(x-xs)/dx;
		jrs = (int)ar;
		if(jrs>nr-2) jrs = nr-2;
		srs = ar-jrs;
		srs0 = 1.0-srs;
		ar = (xg>=x)?(xg-x)/dx:(x-xg)/dx;
		jrg = (int)ar;
		if(jrg>nr-2) jrg = nr-2;
		srg = ar-jrg;
		srg0 = 1.0-srg;
		sigp = ((xs-x)*(xg-x)>0)?1.0:-1.0;
		zpt[ix] = fzt+(nzt-1)*dzt;

		for(iz=izt0; iz<nzt; ++iz){
			sigs = srs0*sigb[jrs][iz]+srs*sigb[jrs+1][iz]; 
			sigg = srg0*sigb[jrg][iz]+srg*sigb[jrg+1][iz]; 
			coss = srs0*cosb[jrs][iz]+srs*cosb[jrs+1][iz]; 
			cosg = srg0*cosb[jrg][iz]+srg*cosb[jrg+1][iz]; 
			ampd = v0*dx*(coss+cosg)*dz;
			if(ampd<0.0) ampd = -ampd;
			if(ls) 
			    ampt[ix][iz] = ampd/sqrt(v0*sigs*sigg);
			else 
			    ampt[ix][iz] = ampd/sqrt(sigs*sigg*(sigs+sigg));

			pd = srs0*pb[jrs][iz]+srs*pb[jrs+1][iz]+sigp 
			     *(srg0*pb[jrg][iz]+srg*pb[jrg+1][iz]);
			if(pd<0.0) pd = -pd;
			if(pd<0.0) pd = -pd;
			temp = 0.5*pd*v0*dx/dz;
			if(temp<1) temp = 1.0;
			if(temp>mzmax) temp = mzmax;
			ampti[ix][iz] = ampt[ix][iz]/(temp*temp);
			zmt[ix][iz] = temp;
			if(pd<pmin && zpt[ix]>fzt+(nzt-1.1)*dzt) 
				zpt[ix] = fzt+iz*dzt;

		}
	}
	
	nxf = (xm-aperx-fx)/dx+0.5;
	if(nxf<0) nxf = 0;
	nxe = (xm+aperx-fx)/dx+0.5;
	if(nxe>=nx) nxe = nx-1;
	
	/* interpolate amplitudes */
	for(ix=nxf; ix<=nxe; ++ix){
		x = fx+ix*dx;
		dis = (xm>=x)?xm-x:x-xm;
		izt0 = (dis*rxz-fzt)/dzt-1;
		if(izt0<0) izt0 = 0;
		if(izt0>=nzt) izt0 = nzt-1;
		iz0 = (dis*rxz-fz)/dz;
		if(iz0<0) iz0 = 0;
		if(iz0>=nz) iz0 = nz-1;

		ax = (x-fxt)/dxt;
		jx = (int)ax;
		ax = ax-jx;
		if(ax<=0.01) ax = 0.;
		if(ax>=0.99) ax = 1.0;
		ax0 = 1.0-ax;
		if(jx>nxte-1) jx = nxte-1;
		if(jx<nxtf) jx = nxtf;

		ar = (xs>=x)?(xs-x)/dx:(x-xs)/dx;
		jrs = (int)ar;
		if(jrs>=nr-1) jrs = nr-2;
		srs = ar-jrs;
		srs0 = 1.0-srs;
		ar = (xg>=x)?(xg-x)/dx:(x-xg)/dx;
		jrg = (int)ar;
		if(jrg>=nr-1) jrg = nr-2;
		srg = ar-jrg;
		srg0 = 1.0-srg;

		for(iz=izt0; iz<nzt; ++iz){
		    tzt[iz] = ax0*tsum[jx][iz]+ax*tsum[jx+1][iz]
				+srs0*tb[jrs][iz]+srs*tb[jrs+1][iz]
				+srg0*tb[jrg][iz]+srg*tb[jrg+1][iz];

		    amp[iz] = ax0*ampt[jx][iz]+ax*ampt[jx+1][iz];
		    ampi[iz] = ax0*ampti[jx][iz]+ax*ampti[jx+1][iz];
		    zm[iz] = ax0*zmt[jx][iz]+ax*zmt[jx+1][iz];
		
		}

		nzp = (ax0*zpt[jx]+ax*zpt[jx+1]-fz)/dz+1.5;
		if(nzp<iz0) nzp = iz0;
		if(nzp>nz) nzp = nz;

		/* interpolate along depth if operater aliasing 	*/
		for(iz=iz0; iz<nzp; ++iz) {
			az = z0+iz*rdz;
			jz = (int)az;
			if(jz>=nzt-1) jz = nzt-2;
			sz = az-jz;
			sz0 = 1.0-sz;
			td = sz0*tzt[jz]+sz*tzt[jz+1];
			at = (td-ft)*odt;
			jt = (int)at;
			if(jt > 0 && jt < nt-1){
			    ampd = sz0*ampi[jz]+sz*ampi[jz+1];
			    mz = (int)(0.5+sz0*zm[jz]+sz*zm[jz+1]);
			    res = at-jt;
			    iz1 = iz+mzmax;
 			    temp = (-migi[ix][iz1-mz]+2.0*migi[ix][iz1]
				-migi[ix][iz1+mz])*ampd;				
			    trace[jt] += (1.0-res)*temp;
			    trace[jt+1] += res*temp;
			}
		}

		/* interpolate along depth if not operater aliasing 	*/
		for(iz=nzp; iz<nz; ++iz) {
			az = z0+iz*rdz;
			jz = (int)az;
			if(jz>=nzt-1) jz = nzt-2;
			sz = az-jz;
			sz0 = 1.0-sz;
			td = sz0*tzt[jz]+sz*tzt[jz+1];
			at = (td-ft)*odt;
			jt = (int)at;
			if(jt > 0 && jt < nt-1){
			    ampd = sz0*amp[jz]+sz*amp[jz+1];
			    res = at-jt;
			    temp = mig[ix][iz]*ampd;
			    trace[jt] += (1.0-res)*temp;
			    trace[jt+1] += res*temp;
			}
		}
	}
					   
		
	
	/* apply half-derivative filter to trace */
	conv(nhd,-lhd,hd,nt,0,trace,nt,0,work1);

	/* convolve wavelet with trace */
 	conv(w->lw,w->iw,w->wv,nt,0,work1,nt,0,trace); 	
	
	/* free workspace */
	free2float(ampt);
	free2float(ampti);
	free1float(tzt);
	free1float(work1);
 	free1float(amp);
 	free1float(ampi);
 	free2float(zmt);
 	free1float(zpt);
 	free1float(zm);
}

